1. Field of the Invention
This invention relates generally to methods for evaluating mathematical functions of many variables, and more specifically to a method for packing and extracting tetrahedra and octahedra within a function domain to approximate the value of a continuous multivariable function.
2. Description of the Prior Art
Many modern-day processes require the evaluation of one or more mathematical functions. These functions frequently involve a multiplicity of variables interrelated by relatively complex mathematical expressions. Oftentimes, the overall efficiency of a process is determined by the relative ease or difficulty with which a function may be evaluated. For instance, in the field of color imaging, it is often desired to convert a color image from a first color space to a second color space. The first color space may be for an image representation stored in a computer file, and the second color space may be for specification of color created on paper by means of printer inks.
During the implementation of a process, it is generally more important to provide an efficient function evaluation technique, as opposed to providing a function evaluation technique optimized for maximum accuracy. For example, during the process of color image conversion, a quick, expedient evaluation of color image functions is much more important than a completely accurate representation of the converted color image. Efficiency is of paramount importance in this context because a common design goal of many image processing systems is to provide rapid processing capability. Accuracy is not a critical design parameter because the human eye cannot distinguish or compensate for minor imaging errors. The characteristics of human visual perception generally allow for the existence of a known bounded error in image rendering. Although the importance of efficiency relative to accuracy has been described in conjunction with color imaging systems, a similar situation exists across a broad spectrum of other process applications as well.
In some process applications, the output of a function may be sampled, but the function itself is unknown. For example, consider the conversion of an image from the color space of a computer file to the surface of a newspaper. The actual colors produced by the printing ink interacting with the paper are difficult or impossible to accurately quantify. The factors resulting in the appearance of a given color on paper are represented by an unknown function, even though the value of the function can be sampled at various intervals. In this scenario, an accurate representation of the imaging information is not at all critical, and emphasis may instead be placed on data processing speed.
Many existing processes involve mathematical functions having extensive domains. It is generally impractical to sample the output of the function at all possible sample values. For example, FIG. 1 illustrates the domain of a mathematical function representing the color of an object in the context of a color image processing system. State of the art color imaging systems represent a color image as a regular array of spots, generally referred to as pixels. Each pixel is assigned a color represented by the coordinates of the color in a three dimensional space. In an additive color environment (i.e., a cathode-ray tube display), any color may be represented by a given combination of the three primary colors of red, green, and blue. With reference to FIG. 1, the value of red may be represented along the X axis 12, green along the Y axis 14, and blue along the Z axis 16. Each coordinate may be of arbitrary precision, but coordinates are generally represented using 8-bit values. In this manner, each pixel may be assigned one of 2.sup.24 different colors. Therefore, the domain of the mathematical function representing pixel color is quite extensive.
The mathematical function accepts input values within the domain of the function, and produces output values corresponding to the input values. The sequence of mathematical operations carried out by the function are determined by the process to be implemented. To print a color image displayed on a cathode-ray tube, the image must be converted from additive to subtractive form. As previously described, the cathode-ray tube image is stored as a pixel-by-pixel representation specifying particular quantities of the colors red, green, and blue. For subtractive color applications, such as printing, the three primary colors are cyan (blue-green), magenta, and yellow with black used as an additional or optional color for increased density ( darkness ) or decreased total ink usage. Accordingly, a function must be determined for the purpose of converting pixel color representations into known quantities of colored printer inks, typically with the amount expressed as an integer in the range of 0 to 255 for each of cyan, magenta, yellow, and black. The function accepts input values for the variables red, green, and blue, and produces output values which represent quantities of cyan, magenta, yellow, and black. Other color spaces in use as either input or output spaces include the colorimetric spaces which represent color based on the tristimulus values that represent a standard observer as defined by the Commission Internationale de l'Eclairage. CIE L*a*b*, CIE L*u*v*, and CIE XYZ are three spaces.
For many processes, it is highly impractical to define a function in analytical form. In the context of color image conversion, it would be very difficult to develop a working analytical model for the purpose of converting a cathode-ray tube image representation into a form suitable for color printing on newspaper media. If the printer inks offered a perfectly linear response, and if the paper was perfectly white, then the function for the conversion of the color image could be specified analytically. However, as a practical matter, the inks have a nonlinear response, and the paper is off-white. Consequently, the print function is best represented by using measurements accumulated for a plurality of print samples. There are 2.sup.24 points in the domain of the print color function which could be sampled, so it would be prohibitively time-consuming and storage-intensive to measure and store all of these values. Rather, the function may be approximated using a smaller set of measured domain values (sample points), and using interpolation to compute the approximate values for all of the other domain values.
FIG. 1 illustrates the concept of taking sample points at regularly-spaced intervals throughout the display-color space, which is the domain of the print-color function. Each sample point is represented by an ordered triplet (R, G, B) which corresponds to values along the red (R), green (G), and blue (B) axes of display-color space. In the present example, samples are taken at points 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46, and 48. Accordingly, the sample points form a rectangular grid throughout the domain of the function, such that the rectangular grid is comprised of a plurality of rectangular volumes 10.
Many existing linear interpolation techniques store the coordinates of the sample points 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, ! 38, 40, 42, 44, 46, and 48 in one or more interpolation tables. The interpolation tables associate each domain sample point coordinate with a function range output value. The domain sample points are generally situated in a first color space, such as red-green-blue, whereas the range output values are situated in a second color space, such as cyan-magenta-yellow.
The accuracy of a given interpolation technique is determined in part by the size of the interpolation tables. In approximating a continuous function, increasing the number of entries in the tables and decreasing the distance between the entries can be used to provide virtually any desired level of function approximation accuracy. However, this approach to increasing system accuracy is very expensive in practice. Table storage and table indexing structures consume large blocks of memory, especially in the context of three-or higher-dimensional functions. Accordingly, it would be desirable to develop an approximating technique which produces accurate answers with relatively small tables.
The overall efficiency of a given approximation technique is often determined by the complexity of mathematical calculations required to perform an approximation. The cost of performing the calculations necessary to produce the result should be minimized. Many prior art systems require the implementation of numerous memory access operations to read table values from a relatively high number of memory locations. One way to minimize the complexity of the mathematical calculations is to reduce the number of locations in the table which the approximation equipment uses to compute the result. By contrast, the complexity of the calculations required to load the interpolation tables is generally not critical. The loading of the interpolation table is carried out infrequently as compared with the number of times the table is used to approximate the function. Therefore, it would be desirable to develop an interpolation technique which minimizes the cost of performing the calculations and the number of memory locations which must be accessed, even if the computational complexity of loading the interpolation table is increased.
State of the art function approximation techniques produce answers which are slightly in error. For system applications such as color image processing, these small errors are generally acceptable. However, some approximation techniques produce discontinuous results in the output range of the function when the domain is presented with continuous input points. Such discontinuous results may produce visible artifacts, because the human eye is more sensitive to differences in color and/or intensity than it is to the absolute value of color and/or intensity. If the magnitude of the discontinuity produced by the approximating scheme exceeds the just-noticeable difference, visible contours will appear in the areas of the image that are supposed to be perfectly smooth. Consequently, it would be desirable to develop an efficient function approximation technique which provides a continuous range output in response to a continuous domain input.
Presently existing function approximation techniques which utilize interpolation tables may be categorized by considering the manner in which the techniques perform each of three subfunctions. The subfunctions are (1) the subdivision of the function domain to select sample domain input: points, (2) the extraction of appropriate sample domain input points corresponding to an arbitrarily selected point from within the function domain, and (3) the actual mathematical interpolation process.
The operation of selecting sample points in the continuous input domain of the function can be termed packing or subdividing. The function is evaluated at these sample points and the evaluation results are saved in a table for subsequent utilization by an interpolation process. The selection of sample points often has a geometrical interpretation: it may be conceptualized as the division of the input domain of the function into a set of contiguous n-dimensional polyhedrons, with the sample points being represented as the polyhedron vertices. For example, a function domain consisting of a three-dimensional rectangular solid can be subdivided into a set of rectangular solids along three mutually orthogonal sets of parallel planes.
Although the division of the function domain into rectangular volumes is a common technique employed by many existing interpolation techniques, it is not an absolute requirement. The function domain may be divided into volumes having nonrectangular geometrical configurations to enhance the efficiency of the interpolation technique in the context of specific system applications. Furthermore, the packing need not be conducted at regular intervals. For example, it may be advantageous to use smaller rectangular solids in areas of the function domain having relatively high amounts of curvature. However, it is necessary that the collection of solids chosen for packing be non-overlapping and completely fill the portion of the function domain for which approximations are desired.
The process of selecting a small number of sample points to be used in computing a given approximation can be termed extraction. Since the sample points typically represent polyhedra vertices, extraction refers to the process of identifying the polyhedron containing the arbitrarily selected point for which an approximation is desired. The arbitrarily selected point may be called the target evaluation point. In the context of an interpolation scheme utilizing tables, the extraction process generally involves extracting the desired function values from a table. The desired function values are the function values at the coordinates of the vertices of the extracted polyhedron.
After the processes of subdivision and extraction have been implemented, the actual interpolation process is conducted. The interpolation operation takes the sample points extracted from the table and uses the function value at these points together with the coordinates of the target evaluation point as inputs for a process which produces the approximate value of the function.
One prior art technique for converting selected points in the domain of a function to values in a function range is known as trilinear interpolation. In the environment of a function, F, having a three-dimensional input domain, the domain is completely filled with rectangular solids. This process is implemented by selecting a series of points along each axis, for example, the x, y, and z axes, as follows: (x.sub.0, x.sub.1, x.sub.2, . . . x.sub.a), (y.sub.0, y.sub.1, y.sub.2, . . . y.sub.b), and (z.sub.0, z.sub.1, z.sub.2, . . . z.sub.c). These series are chosen such that x.sub.i =x.sub.0 +i*(x.sub.a -x.sub.0)/a, y.sub.j =y.sub.0 +j*(y.sub.b -y.sub.0)/b, and z.sub.k =z.sub.0 +k*(z.sub.c -z.sub.0)/c. In this manner, the sample points are (x.sub.i, y.sub.j, z.sub.k) for 0.ltoreq.i.ltoreq.a, 0.ltoreq.j.ltoreq.b, and 0.ltoreq.k.ltoreq.c. Prior to the step of function approximation, the value of the function F(x,y,z) is measured at each of these sample points.
To approximate the function at a target evaluation point (r,s,t) selected from the input domain of the function, the following procedure is used. The function is approximated for the target point (r,s,t) as follows: ##EQU1## where x.sub.i .ltoreq.r.ltoreq.x.sub.i+1, y.sub.j .ltoreq.s.ltoreq.y.sub.j+1, z.sub.k .ltoreq.t.ltoreq.z.sub.k+1
d.sub.r =(r-x.sub.i)/(x.sub.i+1 -x.sub.i) PA1 d.sub.s =(s-y.sub.j)/(y.sub.j+1 -y.sub.j) PA1 d.sub.t =(t-z.sub.k)/(z.sub.k+1 -z.sub.k) PA1 1. The point lies on a straight line connected the two nearest sample points. PA1 2. The point lies on the interior of an equilateral triangle formed by connecting the three nearest sample points. PA1 3. The point lies in the interior of a rectangular tetrahedron formed by connecting the four nearest sample points. PA1 4. The point lies in the interior of a regular octahedron formed from six sample points such that at least three of the points are the nearest sample points and the interior of the regular octahedron does not include any other sample points.
With reference to the formulas set forth above, the evaluation of the function at target evaluation point (r,s,t) requires the use of eight different sample points in the interpolation calculations. Furthermore, at least eight multiplications and seven additions are required per approximation. This interpolation process uses rectangular solids for packing, extracts the rectangular solid that encloses the target point, and uses trilinear interpolation of the eight vertices to calculate the approximation.
Some prior art interpolation systems improve upon the basic method of trilinear interpolation by storing values for d.sub.r, d.sub.s, and d.sub.t, and/or values for various combinations of terms in the above formulas. However, even if some of the formula terms are calculated and stored in a table prior to function approximation, the overall computational efficiency of this trilinear interpolation method is limited by the fact that sample values for eight vertices must be used to calculate each approximation. An important factor in determining the overall effectiveness of linear interpolation methods is the packing technique used to divide the function domain. A packing technique which divides the function domain into geometrical volumes having a minimum number of vertices would greatly simplify the required calculations, as compared to the scheme described above which uses eight vertices per interpolation. What is needed is an improved packing technique for dividing the function domain into geometrical volumes. In this manner, the interpolation process will calculate a function approximation using a minimum number of sample points per calculation.
One exemplary prior art interpolation process which reduces the number of vertices required for a function approximation is described in U.S. Pat. No. 4,477,833 issued to Clark. For purposes of explanation, it is assumed that the function domain is partitioned into unit divisions in each dimension, such that x.sub.i+1 -x.sub.i =1. The domain packing algorithm of Clark uses rectangular solids in a manner substantially identical to that of the standard trilinear method described above. The polyhedral extraction algorithm disclosed in Clark returns four sample points per target evaluation point. The four sample points are a subset of the eight points returned in the standard trilinear method.
The four sample points used in the Clark method are chosen as follows for a target evaluation point (r,s,t). Let (x.sub.i, y.sub.j, z.sub.k) be the sample point such that x.sub.i r.ltoreq.x.sub.i+1, y.sub.j .ltoreq.s.ltoreq.y.sub.j+1, and z.sub.k .ltoreq.t.ltoreq.z.sub.k+1. Define d.sub.r as (r-x.sub.i)/(x.sub.i+1 -x.sub.i), d.sub.s as (s-y.sub.j)/(y.sub.j+1 -y.sub.j), and d.sub.t as (t-z.sub.k)/(z.sub.k+1 -z.sub.k). Then, d.sub.r, d.sub.s, and d.sub.t are used to determine the dominant component, i.e., the maximum d, and the second dominant component.
Define "primary" to be the vertex of the subcube enclosing the target point that corresponds to the dominant component, and primeWeight as follows: primary=(x.sub.i+1, y.sub.j, z.sub.k) and primeWeight=d.sub.r if d.sub.r &gt;d.sub.s and d.sub.r &gt;d.sub.t ; primary=(x.sub.i, y.sub.j+1, z.sub.k) and primeWeight=d.sub.s if d.sub.s &gt;d.sub.r and d.sub.s &gt;d.sub.t ; primary=(x.sub.i, y.sub.j, z.sub.k+1) and primeWeight=d.sub.t if d.sub.t &gt;d.sub.r and d.sub.t &gt;d.sub.s.
Define "secondary" as the corner adjacent to the primary and in the direction of the second dominant component. Consider the definition for primary=(x.sub.i+1, y.sub.j, z.sub.k) where d.sub.r is the dominant component. Secondary=(x.sub.i+1, y.sub.j+1, z.sub.k) secondWeight=d.sub.s, and lastWeight=d.sub.t if d.sub.s &gt;d.sub.t ; Secondary=(x.sub.i+1, y.sub.j, z.sub.k+1) secondWeight=d.sub.t, and lastWeight=d.sub.s if d.sub.t &gt;d.sub.s.
In accordance with the Clark patent, the approximation for the function value is calculated as follows: ##EQU2## In this manner, the Clark process uses rectangular solids for packing the function domain. An extraction algorithm extracts four sample points, and a directional interpolation algorithm of the four sample points approximates the function value corresponding to a target evaluation point.
The Clark process yields approximation results which contain a certain amount of error. The extent of the error is partially dependent upon the particular extraction method employed, the nature of the function to be approximated, and the domain of interest. Therefore, in some applications, the Clark method may provide acceptable results, whereas, when applied to other applications, the Clark method may prove wholly unsatisfactory. What is needed is a function approximation method which utilizes the relatively low number of sample points employed by the Clark method, yet offers enhanced accuracy in applications where the Clark method would not be well-suited.
Another exemplary prior art is U.S. Pat. No. 4,275,413 issued to Sakamoto, et al, which describes a tetrahedral volumetric interpolation method for color correction. This method uses a domain packing algorithm similar to Clark. The extraction process returns 4 points, each of which is a vertex of the rectangular solid, or an average of 4 or 8 of the vertices of the rectangular solid. Their averaged points are the centers of a face of the rectangular solid (i.e., the average of 4 corners of a rectangle) or the center of the rectangular solid (i.e., the average of the 8 vertices of the rectangular solid).